---
title: "Analysis of Princeton Trade Study Results"
author: "Robert Kubinec"
output:
  html_document:
    df_print: paged
---

```{r setup, include=FALSE}
require(tidyr)
require(dplyr)
require(ggplot2)
require(readr)
require(lubridate)
require(stringr)
require(margins)
require(qualtRics)
require(ggthemes)
require(googlesheets4)
require(modelsummary)
require(RPostgres)

# load data

# this is for the raw data, data in repo is stripped of any identifying information

# qualtrics <- read_survey("data/U.S.+Firms+Outcome+Survey+-+API+Based_October+5,+2020_08.17.csv") %>% mutate(survey_type="Facebook") %>% 
#   filter(Status!="Survey Preview")
# 
# ks <- read_survey("data/U.S.+Firms+Outcome+Survey+-+KC+Chamber_October+5,+2020_06.18.zip") %>% mutate(survey_type="Kansas City Chamber") %>% 
#   filter(Status!="Survey Preview")
# 
# # note this doesn't contain any responses
# 
# phone <- read_survey("data/U.S.+Firms+Outcome+Survey+-+Phone_October+5,+2020_08.24.csv") %>% mutate(survey_type="Orbis Phone") %>% 
#   filter(Status!="Survey Preview")
# 
# email <- read_survey("data/U.S.+Firms+Outcome+Survey_October+5,+2020_08.26.csv") %>% mutate(survey_type="Orbis Email") %>% 
#   filter(Status!="Survey Preview")
# 
# # need facebook valid variable
# 
# valid_fb <- googlesheets4::read_sheet("https://docs.google.com/spreadsheets/d/1SlyYtROA6K1HIpx9MQ1QRhGyDpIkwfwebKGHguZ3b98/edit#gid=1263452815",
#                                       sheet="Facebook")
# 
# qualtrics <- left_join(qualtrics,select(valid_fb, ResponseId,Valid))
# 
# # only those who Lindsay paid
# valid_zero <- filter(qualtrics, Valid==0)
# # additinal set that Lindsay padi
# valid_one <- filter(qualtrics,Valid==1)
# # add in those considered valid
# valid_maybe <- filter(qualtrics,Valid==0.5)
# 
# valid_most <- filter(qualtrics,Valid>0)
# 
# # new survey responses
# 
# new_survey <- filter(qualtrics,StartDate>ymd("2020-10-01"))

# bind together into one dataset

# combined_data <- bind_rows(list(`Full Dataset`=qualtrics,
#                                 `Not Valid`=valid_zero,
#                                 `Valid`=valid_one,
#                                 `Maybe Valid`=valid_maybe,
#                                 `Excluding Invalid`=valid_most,
#                                 `New Survey`=new_survey),.id = "Type")

# select(combined_data, -IPAddress, -LocationLatitude, -LocationLongitude, -email_7, -email_8,
#        -email, -matches("invite")) %>%
#   saveRDS(file="data/combined_anon.rds")

combined_data <- readRDS("data/combined_anon.rds")

# combine paid valid all with ks survey, orbis survey and email data
# for final analysis



# to_analyze <- bind_rows(valid_most,ks,email,new_survey) %>% 
#   filter(Status!="Survey Preview")

# select(to_analyze, -IPAddress, -LocationLatitude, -LocationLongitude, -email_7, -email_8,
#        -email, -matches("invite")) %>% 
#   saveRDS(file="data/to_analyze_anon.rds")

# load anonymized datasets

to_analyze <- readRDS("data/to_analyze_anon.rds")


knitr::opts_chunk$set(warning=FALSE,message=FALSE)

```

This file contains the code to reproduce analyses in "Tariff and Corporate Political Activity: A Survey on U.S. Businesses", forthcoming in *Business and Politics.* For any questions about the code, please contact Robert Kubinec at bobkubinec@gmail.com.

Please note that this script produces additional analyses that are not included in the published paper.


```{r getnaics,include=F}

# BEA industry data

beapretties <- readRDS("beapretties.rds")

```

```{r getnaics2,include=F}

# see if we can join

analyze_tariffs <- left_join(to_analyze,beapretties, by=c("naics_2"="naicsdesc")) %>% 
  group_by(ResponseId) %>% 
  mutate(number=coalesce(number,0),
         prod_sum=sum(number),
         total=sum(number>0),
         mean_rate=mean(rate[number>0]),
         mean_prop=mean(tariff_prop_weighted[number>0]),
        prod_sum=prod_sum/100,
        mean_rate=mean_rate/100) %>% 
  distinct(ResponseId,prod_sum,mean_rate,mean_prop,.keep_all = T) %>% 
  ungroup %>% 
  ungroup %>% 
  mutate(prod_sum_orig=prod_sum,
         prod_sum=ifelse(condition %in% c("dynamic_static","static") & prod_sum==0,
                        (-1*min(prod_sum[prod_sum>0],na.rm=T)),
                        prod_sum))

# plot treatment distribution

analyze_tariffs %>% 
  mutate(Treatment=(condition!="control_email"),
         Treatment=factor(Treatment, labels=c("Control","Treatment"))) %>% 
  filter(!is.na(condition)) %>% 
  ggplot(aes(x=prod_sum*100)) +
  geom_histogram() +
  theme_tufte() +
  facet_wrap(~Treatment,ncol=1,scales="free_y") +
  labs(x="Number of Products with Tariffs in Respondent's Industry",
       y="Count of Respondents")

ggsave("treathist.png",width=5, height=3)

# table of respondents by tariff category

```



We have different datasets depending on the level of validation we used. The validation column is unfortunately raw text so right now I am only using a simple text match in the `badresponse` column for the strings `bot` and `fake` and also including those who were paid in the first round. I also have a dataset of only those who were paid (first round) and all who were paid in the first and second rounds. Finally, I am also performing analyses for the full dataset as a reference point.

# Descriptive Analysis

First we can do some descriptive results. First for information about whether or not respondents know about tariffs by dataset type:

```{r know_all}
combined_data %>% 
  ggplot(aes(x=factor(know_trade_1,
                      levels=c("0","1","2","3","4","5","6","7","8","9","10")))) +
  geom_bar() +
  ggtitle("On a scale of 0 (very low) to 10 (very high),\nhow much information do you already have about how your \ncompany has been affected by the trade war?") +
  facet_wrap(~Type,scales="free_y") +
  xlab("")

ggsave("compare_info.png")

```


It is not clear if there are any patterns here across the respondent validity. There are distinct modes at 0, 7, and 10. The full dataset has a higher mode at 7 than the other datasets.

We can see a similar plot for the trade war's effect on the respondent's company:

```{r trade_war_effect}
combined_data %>% 
  ggplot(aes(x=factor(hurt_trade_1,
                      levels=c("0","1","2","3","4","5","6","7","8","9","10")))) +
  geom_bar() +
  ggtitle("On a scale of 0 (harmed) to 10 (helped), how has your company been affected by the trade tariffs that went into effect in over the past year?") +
  facet_wrap(~Type,scales="free_y") +
  xlab("")

ggsave("compare_trade_war_hurt.png")
```

Here the distributions across the categories are even more equal regardless of respondent validity. There is a distinct mode at 5 and a lesser mode at 10 (people who are very positive about the trade war). 

Here are the same descriptive statistics for several other pre-treatment questions:

```{r other_stats}
combined_data %>% 
  ggplot(aes(x=factor(effective,
                      levels=c("Not effective at all",
                               "Somewhat ineffective",
                               "Neither effective nor ineffective",
                               "Somewhat effective",
                               "Very effective")))) +
  geom_bar() +
  ggtitle("On the whole, would your company consider\nit an effective strategy to work\nwith other companies to advocate for a common trade policy?") +
  facet_wrap(~Type,scales="free_y") +
  xlab("") + 
  theme(axis.text.x = element_text(angle=90))

combined_data %>% 
  ggplot(aes(x=pol_action)) +
  geom_bar() +
  ggtitle("In the last five years, has your company ever taken\npolitical action on any issue?") +
  facet_wrap(~Type,scales="free_y") +
  xlab("")

combined_data %>% 
  ggplot(aes(x=assoc)) +
  geom_bar() +
  ggtitle("Does your company participate\nin trade and business associations that advocate\nfor your industry?") +
  facet_wrap(~Type,scales="free_y") +
  xlab("")

# bind_rows(list(Kansas=select(ks,pol_culture_1),
#                Facebook=select(valid_most,pol_culture_1)),.id="Type") %>% 
#   filter(!is.na(pol_culture_1)) %>% 
#   ggplot(aes(x=factor(pol_culture_1,
#                       levels=c("Very Liberal",
#                                "Liberal",
#                                "Moderate",
#                                "Conservative",
#                                "Very Conservative")))) +
#   geom_bar() +
#   ggtitle("How would you evaluate your company’s\npolitical culture among managers?") +
#   facet_wrap(~Type,scales="free_y") +
#   xlab("") + 
#   ylab("No. Respondents") +
#   theme_tufte() +
#   theme(axis.text.x = element_text(angle=90))
# 
# ggsave("pol_culture_managers_compare.png")

combined_data %>% 
  ggplot(aes(x=factor(pol_culture_2,
                      levels=c("Very Liberal",
                               "Liberal",
                               "Moderate",
                               "Conservative",
                               "Very Conservative")))) +
  geom_bar() +
  ggtitle("How would you evaluate your company’s\npolitical culture among rank and file employees?") +
  facet_wrap(~Type,scales="free_y") +
  xlab("") + 
  theme(axis.text.x = element_text(angle=90))

```

Again, we're not seeing very strong effects due to validation type. The full dataset (absent bots) are more likely to be politically active. 

# More Descriptive Stuff

In this plot, we look at the total numbers of products with tariffs by different industries.

```{r tarifftable}

beapretties %>% 
  filter(!is.na(number)) %>% 
  distinct(industry_3dig_desc,number,tariff) %>% 
  mutate(industry_3dig_desc=str_trunc(industry_3dig_desc,30,side="right"),
         tariff=case_match(tariff,
                           "Washing Machine"~"Washing Machine (25%)",
                           "Aluminum"~"Aluminum (10%)",
                           "Steel"~"Steel 25%",
                           "Section 301 (China)"~"Section 301 (17.5%)")) %>% 
  ggplot(aes(y=reorder(industry_3dig_desc,number),x=number)) +
  geom_col(aes(fill=tariff)) + 
  scale_fill_grey(name="Tariff",guide=guide_legend(nrow=2)) + 
  ggthemes::theme_clean() + 
  theme(legend.position = "bottom",
        legend.text = element_text(size=10),
        legend.title.position = "top") + 
  scale_y_discrete(guide=guide_axis(check.overlap = T)) +
  labs(y="",x="Count of Input Products with Tariffs")

ggsave("tariff_count_plot.pdf",width=6,height=7)
 
```
Here we can look at the reasons that people did not select any political outcomes for the trade war.

```{r barnone}

to_analyze %>% 
  filter(!is.na(none_reason),!is.na(outcomtype)) %>% 
  mutate(none_reason=ifelse(none_reason=="Even if the trade war is harmful in the short run, we think in the long run it will benefit the U.S. economy.","Even if the trade war is harmful in the short run, we think it will benefit the U.S. economy in the long run.",none_reason),
         none_reason=stringr::str_wrap(none_reason, width=50),
         none_reason=factor(none_reason, levels = c("Something else:",setdiff(unique(none_reason), "Something else:"))),
         outcomtype=str_wrap(outcomtype, width=15)) %>% 
ggplot(aes(x = none_reason)) +
  geom_bar() +
  facet_wrap(~ outcomtype,nrow=1,scales="free_x") +
  labs(x = "",
       y = "Count of Responses") +
  ggthemes::theme_clean() +
  theme(axis.text.x = element_text(angle = 40, hjust = 1)) +
  coord_flip()

ggsave("reason_war.pdf",width=6,height=6)

```


# Treatment

Now I will look at the treatment. In our current survey we have treatments for taking supporting and opposing stances in response to the data shown. We can count the number of possible actions taken and use the treatment condition as a predictor to see what the results are.

After munging the data, I first run a model with the treatments predicting opposing the trade war as a binary variable where 1 equals the respondent selecting at least one opposition outcome. Whether or not the respondent was in the paid/and or not bot category is included as a control.

```{r rework_treatment}
# need to split the treatment variable
to_analyze <- to_analyze %>% 
  mutate(support_fb=grepl("Facebook",outcome_support),
         support_congress_donate=grepl("Donate to Congress",outcome_support),
         support_petition=grepl("Petition",outcome_support),
         support_invite=grepl("Invite",outcome_support),
         oppose_fb=grepl("Facebook",outcome_oppose),
        oppose_congress_donate=grepl("Donate to Congress",outcome_oppose),
        oppose_petition=grepl("Petition",outcome_oppose),
        oppose_invite=grepl("Invite",outcome_oppose),
        oppose_congress_ask=grepl("Ask your Congress",outcome_oppose),
        oppose_gov=grepl("governor",outcome_oppose),
        support_any=as.numeric(support_fb|support_congress_donate|support_petition),
        oppose_any=as.numeric(oppose_fb|oppose_congress_donate|oppose_congress_ask|oppose_petition|oppose_congress_ask|oppose_gov),
        oppose_any_view=as.numeric(outcomtype=="Oppose")) %>% 
  filter(!is.na(manager))

supp_mod_list <- list(Facebook=glm(support_fb~condition + manager,family="binomial",data=to_analyze),
                 Congress=glm(support_congress_donate~condition,family="binomial",data=to_analyze),
                 Petition=glm(support_petition~condition,family="binomial",data=to_analyze),
                 Invite=glm(support_invite~condition,family="binomial",data=to_analyze),
                 Any=glm(support_any~condition,family="binomial",data=to_analyze))

opp_mod_list <- list(Facebook=glm(oppose_fb~condition,family="binomial",data=to_analyze),
                 Congress=glm(oppose_congress_ask~condition,family="binomial",data=to_analyze),
                 Petition=glm(oppose_petition~condition,family="binomial",data=to_analyze),
                 Invite=glm(oppose_invite~condition,family="binomial",data=to_analyze),
                 Governor=glm(oppose_gov~condition,family="binomial",data=to_analyze),
                 Any=glm(oppose_any~condition,family="binomial",data=to_analyze))

```

Let's export a table of how many people responded to the prompt in total:

```{r}

library(tinytable)

# frequencies of treatment/control

to_analyze %>% 
  gather(key="Outcome",value="Response",matches("oppose\\_"),-action_tariff) %>% 
  filter(Outcome !="oppose_any_view") %>% 
  group_by(action_tariff,Outcome, Response) %>% 
  count %>% 
  group_by(Outcome) %>% 
  mutate(Response=factor(Response, labels=c("Not Selected","Selected"),
                         levels=c(0,1))) %>% 
  ungroup %>% 
  mutate(Outcome=case_match(Outcome,
                            "oppose_any"~"Any Selection",
                            "oppose_congress_ask"~"Contact Congress",
                            "oppose_congress_donate"~"Donate to Congress",
                            "oppose_fb"~"Join Facebook Group",
                            "oppose_gov"~"Donate to Governor",
                            "oppose_petition"~"Sign Petition",
                            "oppose_invite"~"Invite Others")) %>% 
    group_by(Outcome,action_tariff) %>% 
  mutate(n=paste0(round(n/sum(n)*100,1),"%")) %>% 
    arrange(Outcome,action_tariff) %>% 
  filter(Response!="Not Selected") %>% 
  select(Outcome, `Prior Tariff Action`="action_tariff",
          `Proportion Choosing Outcome`="n") %>% 
  tt(caption="Proportion Selecting Experiment Outcomes by Prior Tariff Action\\label{tbl-outcome-prop}",
    notes = "Table shows proportion selecting different kinds of experimental outcomes given whether the respondents had said they had previously taken action about tariffs.",
    width=.8) %>% 
  format_tt(escape=TRUE) %>% 
  save_tt("exp_prop.tex",overwrite = T)


```

## Reviewer 2

Now redo the same models but collapse treatments:

```{r collapse}

to_analyze <- mutate(to_analyze,condition_collapse=ifelse(condition %in% c("static","dynamic",
                                                                  "dynamic_static"),
                                                 "Treatment","Control"),
                     condition_collapse=factor(condition_collapse,levels=c("Control","Treatment")))

supp_mod_list_collapse <- list(Facebook=glm(support_fb~condition_collapse,family="binomial",data=to_analyze),
                 Congress=glm(support_congress_donate~condition_collapse,family="binomial",data=to_analyze),
                 Petition=glm(support_petition~condition_collapse,family="binomial",data=to_analyze),
                 Invite=glm(support_invite~condition_collapse,family="binomial",data=to_analyze),
                 Any=glm(support_any~condition_collapse,family="binomial",data=to_analyze))

opp_mod_list_collapse <- list(Facebook=glm(oppose_fb~condition_collapse,family="binomial",data=to_analyze),
                 Congress=glm(oppose_congress_ask~condition_collapse,family="binomial",data=to_analyze),
                 Petition=glm(oppose_petition~condition_collapse,family="binomial",data=to_analyze),
                 Invite=glm(oppose_invite~condition_collapse,family="binomial",data=to_analyze),
                 Governor=glm(oppose_gov~condition_collapse,family="binomial",data=to_analyze),
                 Any=glm(oppose_any~condition_collapse,family="binomial",data=to_analyze))

```


Need to do an interaction between prior tariff involvement & the treatment.

```{r prior_tariff_int1}

opp_mod_tarrifp_list <- list(Facebook=glm(oppose_fb~condition*action_tariff,family="binomial",data=to_analyze),
                 Congress=glm(oppose_congress_ask~condition*action_tariff,family="binomial",data=to_analyze),
                 Petition=glm(oppose_petition~condition*action_tariff,family="binomial",data=to_analyze),
                 Invite=glm(oppose_invite~condition*action_tariff,family="binomial",data=to_analyze),
                 Governor=glm(oppose_gov~condition*action_tariff,family="binomial",data=to_analyze),
                 Any=glm(oppose_any~condition*action_tariff,family="binomial",data=to_analyze))


modelsummary(opp_mod_tarrifp_list,stars = TRUE,
             coef_rename = c(conditiondynamic="Dynamic",
                             conditiondynamic_static="Dynamic + Static",
                             conditionstatic="Static",managerYes="Manager",
                             condition_collapseTreatment="Collapsed",
                             action_tariffYes="Tariff Action"),output = "tinytable",
             digts=3,notes="Each model is one of the pre-registered outcomes in opposition to the trade war. Tariff Action is whether the respondent's firm had previously taken action on tariffs.",
     caption="Interaction of Disaggregated Treatments with Prior Tariff Action",width=1) %>% 
  save_tt("tariff_action_disag.tex",overwrite=TRUE)

opp_mod_tarrifp_collapse_list <- list(Facebook=glm(oppose_fb~condition_collapse*action_tariff,family="binomial",data=to_analyze),
                 Congress=glm(oppose_congress_ask~condition_collapse*action_tariff,family="binomial",data=to_analyze),
                 Petition=glm(oppose_petition~condition_collapse*action_tariff,family="binomial",data=to_analyze),
                 Invite=glm(oppose_invite~condition_collapse*action_tariff,family="binomial",data=to_analyze),
                 Governor=glm(oppose_gov~condition_collapse*action_tariff,family="binomial",data=to_analyze),
                 Any=glm(oppose_any~condition_collapse*action_tariff,family="binomial",data=to_analyze))


modelsummary(opp_mod_tarrifp_collapse_list,stars = TRUE,
             coef_rename = c(conditiondynamic="Dynamic",
                             conditiondynamic_static="Dynamic + Static",
                             conditionstatic="Static",managerYes="Manager",
                             condition_collapseTreatment="Collapsed",
                             action_tariffYes="Tariff Action"),output = "tinytable",
             digts=3,notes="Each model is one of the pre-registered outcomes in opposition to the trade war. Tariff Action is whether the respondent's firm had previously taken action on tariffs.",
     caption="Interaction of Disaggregated Treatments with Prior Tariff Action",width=1) %>% 
  save_tt("tariff_action_collapse.tex",overwrite = T)


```


Apparently we need a model in which we "control" for pre-treatment covariates.

```{r prior_tariff_int2}

opp_mod_tarrifp2_list <- list(Facebook=glm(oppose_fb~condition + action_tariff,family="binomial",data=to_analyze),
                 Congress=glm(oppose_congress_ask~condition + action_tariff,family="binomial",data=to_analyze),
                 Petition=glm(oppose_petition~condition + action_tariff,family="binomial",data=to_analyze),
                 Invite=glm(oppose_invite~condition + action_tariff,family="binomial",data=to_analyze),
                 Governor=glm(oppose_gov~condition + action_tariff,family="binomial",data=to_analyze),
                 Any=glm(oppose_any~condition + action_tariff,family="binomial",data=to_analyze))


modelsummary(opp_mod_tarrifp2_list,stars = TRUE,
             coef_rename = c(conditiondynamic="Dynamic",
                             conditiondynamic_static="Dynamic + Static",
                             conditionstatic="Static",managerYes="Manager",
                             condition_collapseTreatment="Collapsed",
                             action_tariffYes="Tariff Action"),output = "tinytable",
             digits=3,notes="Each model is one of the pre-registered outcomes in opposition to the trade war. Tariff Action is whether the respondent's firm had previously taken action on tariffs.",
     caption="Disaggregated Treatments with Prior Tariff Action as a Control",width=1) %>% 
  save_tt("tariff_action_control.tex",overwrite=TRUE)


```


# Extra Models That Are Only in the Appendix


```{r allsupp}

modelsummary(supp_mod_list,output="all_supp.tex",stars=T,
             fmt=3,
             title="Disaggregated Treatments for Any Support Trade War Item Selected",
             coef_rename = c(conditiondynamic="Dynamic",
                             conditiondynamic_static="Dynamic + Static",
                             conditionstatic="Static",managerYes="Manager",
                             condition_collapseTreatment="Collapsed"))

```

```{r allopp}

modelsummary(opp_mod_list,output="all_opp.tex",stars=T,
             fmt=3,
             title="Disaggregated Treatments for Any Oppose Trade War Item Selected",
             coef_rename = c(conditiondynamic="Dynamic",
                             conditiondynamic_static="Dynamic + Static",
                             conditionstatic="Static",managerYes="Manager",
                             condition_collapseTreatment="Collapsed"))

```


```{r collapseopp}

modelsummary(opp_mod_list_collapse,output="collapse_opp.tex",stars=T,
             title="Collapsed Treatments for Any Oppose Trade War Item Selected",
             fmt=3,
             coef_rename = c(conditiondynamic="Dynamic",
                             conditiondynamic_static="Dynamic + Static",
                             conditionstatic="Static",managerYes="Manager",
                             condition_collapseTreatment="Collapsed"))

```

```{r collapsesupp}

modelsummary(supp_mod_list_collapse,output="collapse_supp.tex",stars=T,
             title="Collapsed Treatments for Any Support Trade War Item Selected",
             fmt=3,
             coef_rename = c(conditiondynamic="Dynamic",
                             conditiondynamic_static="Dynamic + Static",
                             conditionstatic="Static",managerYes="Manager",
                             condition_collapseTreatment="Collapsed"))

```



Also collapse by counting outcomes:

```{r count}

do_count <- function(outcome) {
  
  outcome_split <- stringr::str_split(outcome,pattern=",")
  
  outcome_split <- sapply(outcome_split, function(l) {
    
   if(is.na(l[[1]])) {
      return(0)
    } else if(l[[1]]=="None of the above") {
      return(0)
    }  else {
      return(length(unique(l)))
    }
  })
  
  outcome_split
  
}

to_analyze <- mutate(to_analyze,count_support=do_count(outcome_support),
                     count_oppose=do_count(outcome_oppose))

supp_mod_count <- list(Oppose=glm(count_oppose~condition,family="gaussian",data=to_analyze),
                               Support=glm(count_support~condition,family="gaussian",data=to_analyze),
                 `Collapse Oppose`=glm(count_oppose~condition_collapse,family="gaussian",data=to_analyze),
                `Collapse Support`=glm(count_support~condition_collapse,family="gaussian",data=to_analyze))


modelsummary(supp_mod_count,stars=T,
             output="count_mods.tex",
             title="Models of Count of Outcome Items Selected",
             coef_rename = c(conditiondynamic="Dynamic",
                             conditiondynamic_static="Dynamic + Static",
                             conditionstatic="Static",managerYes="Manager",
                             condition_collapseTreatment="Collapsed"))

```





# Interaction with Pre-treatment Beliefs

Crucially, we wanted to test whether or not respondents' pre-existing beliefs matter for the effect of the treatment. We in particular hypothesized a curvilinear relationship between beliefs and the outcome. Here I will only analyze the `to_analyze` data.

```{r treatment_int_hurt}

# rework paid-valid-all data

to_analyze <- to_analyze %>%   
    mutate(hurt_trade_1=as.numeric(hurt_trade_1))

# run linear interaction model

linear_mod <-  glm(oppose_any~condition_collapse*hurt_trade_1,data=to_analyze)

# run curvilinear

out_mod_curv <- glm(formula=oppose_any ~ condition_collapse*hurt_trade_1 + condition_collapse*(I(hurt_trade_1^2)),
      data=to_analyze,family="binomial") 

# look at marginal effects

margin_curv_oppose <- margins(out_mod_curv,at = list(hurt_trade_1=seq(0,10,by=1))) %>%
    summary

margin_curv_oppose %>% 
  filter(factor=="condition_collapseTreatment") %>% 
  ggplot(aes(y=AME,x=hurt_trade_1)) +
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="blue",alpha=0.5) +
  geom_line(linetype=2,aes(group=factor)) +
  geom_hline(linetype=3,yintercept=0) +
  theme_tufte() +
  xlab(stringr::str_wrap("On a scale of 0 (harmed) to 10 (helped), how has your company been affected by the trade tariffs that went into effect in over the past year?"))

ggsave("quad_oppose.png")

# same but for support

# run linear interaction model

linear_mod_supp <-  glm(support_any~condition_collapse*hurt_trade_1,data=to_analyze)

# run curvilinear

out_mod_curv_supp <- glm(formula=support_any ~ condition_collapse*hurt_trade_1 + condition_collapse*(I(hurt_trade_1^2)),
      data=to_analyze,family="binomial") 

# look at marginal effects

margin_curv_support <- margins(out_mod_curv_supp,at = list(hurt_trade_1=seq(0,10,by=1))) %>%
    summary

margin_curv_support %>% 
  filter(factor=="condition_collapseTreatment") %>% 
  ggplot(aes(y=AME,x=hurt_trade_1)) +
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="blue",alpha=0.5) +
  geom_line(linetype=2,aes(group=factor)) +
  geom_hline(linetype=3,yintercept=0) +
  theme_tufte() +
  xlab(stringr::str_wrap("On a scale of 0 (harmed) to 10 (helped), how has your company been affected by the trade tariffs that went into effect in over the past year?"))

ggsave("quad_support.png")

```

Interestingly, there are non-linear patterns, though not necessarily those we expect. For the static condition, those who believe that the trade war is very helpful or very hurtful are more likely to oppose the trade war conditional on receiving the treatment:

```{r margins_plot}
# cycle over values of the predictors calculated over average predictions for the data

over_beliefs_static <- prediction(out_mod_curv,at=list(condition_collapse=c("Control","Treatment"),
                                          hurt_trade_1=seq(0,10,by=0.1)))
sum_pred <- summary(over_beliefs_static)

sum_pred %>% 
  ggplot(aes(y=Prediction,x=`at(hurt_trade_1)`)) +
  geom_line(aes(colour=`at(condition_collapse)`,linetype=`at(condition_collapse)`)) +
  geom_ribbon(aes(ymin=lower,ymax=upper,fill=`at(condition_collapse)`),alpha=0.5) +
  #geom_ribbon(data=just_cont,aes(ymin=lower,ymax=upper),fill="red",alpha=0.5) +
  geom_vline(xintercept = 5,linetype=3) +
  #facet_wrap(~`at(condition_collapse)`,ncol=1) +
  theme_minimal() +
  xlab("Scale of 1 to 10 for Harmed by Trade (10 the Most Helped)") + 
  ylab("Probability of Opposing Trade War") + 
  scale_y_continuous(labels=scales::percent) + 
  guides(fill=guide_legend(title=""),
         linetype=guide_legend(title=""),
         colour=guide_legend(title="")) +
  theme(panel.grid=element_blank()) +
  ggtitle("Effect of Treatment on Probability of Selecting Outcome Opposing Trade War",
          subtitle = "Subset by Question: How Has Your Company Been Affected\nby the Trade Tariffs That Went Into Effect Last Year?")

ggsave("treat_hurt_trade.png")
```


We can also investigate an interaction concerning *how much people know about the trade war*. 


```{r know_int}

# run curvilinear
to_analyze$know_trade_1 <- as.numeric(to_analyze$know_trade_1)

out_mod_know <- glm(formula=oppose_any ~ condition_collapse*know_trade_1*hurt_trade_1,
      data=to_analyze,family="binomial") 

```

```{r margins_biv}
# cycle over values of the predictors calculated over average predictions for the data

all_trade <- expand.grid(know_trade_1=seq(0,10,by=0.1),
                         hurt_trade_1=seq(0,10,by=0.1))

over_beliefs_biv1 <- parallel::mclapply(1:nrow(all_trade), function(s) {
  
  print(s)
  
  out_d <- summary(prediction(out_mod_know,
                                  at=list(condition_collapse=c("Control","Treatment"),
                                          know_trade_1=all_trade$know_trade_1[s],
                                          hurt_trade_1=all_trade$hurt_trade_1[s])))
  
  return(out_d)
  
},mc.cores = 3) 
  

sum_pred <- bind_rows(over_beliefs_biv1)

# just_cont <- filter(sum_pred,`at(condition)`=="control_email") %>% 
#   select(-`at(condition)`)

filter(sum_pred,`at(condition_collapse)`!="control_email") %>% 
  mutate(condition_collapse=recode(`at(condition_collapse)`,
                          treatment="Treatment",
                          control="Control")) %>% 
  ggplot(aes(y=`at(know_trade_1)`,x=`at(hurt_trade_1)`)) +
  geom_raster(aes(fill=Prediction)) +
  scale_fill_distiller(type='div') +
  facet_wrap(~condition_collapse) +
  theme_tufte() +
  xlab("Helped by Trade War") + 
  ylab("Knowledge of Trade War") + 
  guides(fill=guide_colorbar(title="Pr(Oppose)")) +
  theme(panel.grid=element_blank())

ggsave("both_know_hurt_trade.png")
```


These results suggest that the more people know about the trade war, the less likely they are to respond to the treatment (which makes sense as the information is less new to them).

Finally we can look at a three-way interaction between the treatment, knowing about the trade war and believing the trade war is harmful. Here I recode the treatment into a binary variable to increase power:


This one is tricky of course to interpret but the sign on the last interaction term suggests that as knowledge about tariffs and beliefs in the beneficial aspects of tariffs both increase, respondents are *more* likely to oppose the tariffs conditional on receiving information. So this result (though imprecise) suggests some initial support for our hypothesis.

# Final Interaction: Political Beliefs

Political beliefs seem to have a strong effect on beliefs about trade.

```{r polmod}

library(marginaleffects)

to_analyze <- mutate(to_analyze,pol_culture_num=as.numeric(factor(pol_culture_1,
                                                     levels=c("Very Conservative",
                                                              "Conservative",
                                                              "Moderate",
                                                              "Liberal",
                                                              "Very Liberal"))),
                     moderate=as.numeric(pol_culture_1=="Moderate"))

pol_managers <- glm(oppose_any ~ condition_collapse*pol_culture_1,data=to_analyze,
                    family="binomial")

pol_rank_file <- glm(oppose_any ~ condition_collapse*pol_culture_2,data=to_analyze,
                    family="binomial")

# get predictions

predict_man <- summary(prediction(pol_managers,at=list(pol_culture_1=c("Very Conservative",
                                                               "Conservative",
                                                               "Moderate",
                                                               "Liberal",
                                                               "Very Liberal"),condition_collapse=c("Control","Treatment")))) %>% 
  mutate(type="Managers") %>% 
  select(estimate=`at(pol_culture_1)`,treatment=`at(condition_collapse)`,everything())

predict_rank <- summary(prediction(pol_rank_file,at=list(pol_culture_2=c("Very Conservative",
                                                               "Conservative",
                                                               "Moderate",
                                                               "Liberal",
                                                               "Very Liberal"),condition_collapse=c("Control","Treatment")))) %>% 
  mutate(type="Rank and File") %>% 
  select(estimate=`at(pol_culture_2)`,treatment=`at(condition_collapse)`,
         everything())

bind_rows(predict_rank,predict_man) %>% 
  ggplot(aes(y=Prediction,x=estimate)) +
  geom_pointrange(aes(ymin=lower,ymax=upper,colour=treatment),position = position_dodge(width=.5)) +
  scale_y_continuous(labels=scales::percent_format(accuracy=1)) +
  theme_tufte() +
  scale_color_colorblind(name="") +
  facet_wrap(~type) +
  scale_x_discrete(guide = guide_axis(n.dodge = 3)) +
  guides(shape=guide_legend(title="")) +
  labs(y="Percent Choosing Oppose Trade War",
       x="")

ggsave("politics_trade_control_treat.png",width=5,height=5)


```


